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SHOCK FITTING APPLIED TO THE PREDICTION OF HIGH-SPEED ROTOR NOISE 


Major John W. Rutherford 

U.S. Array Aviation Research and Technology Activity 
Aeroflightdynamics Directorate 
Ames Research Center 
Moffett Field, California 94035, U.S. A. 


ABSTRACT 


A shock fitting method applied to the transonic small disturbance 
(TSD) potential equation is described. This method is then applied to a 
simple, two-dimensional (2-D) rotating disturbance which is analogous to 
a shock radiating from the tip of a rotor blade in high-speed hover. A 
comparison is made between the results of this method and the more 
standard shock-capturing method. This comparison makes it clear that 
the effect of the results on the acoustic signature of the 2-D model is 
significant, and similar results can be expected when the method is 
extended to the three-dimensional (3-D) case. 


1 . NOMENCLATURE 


a speed of sound, m/sec 

c bump chord length, m 

R radius of model cylinder, m 

r cylindrical radial distance, m 

t time, sec 

y ratio of specific heats 

0 angle measured in cylindrical coordinates 

4> velocity potential, m 2 /sec 

ft angular velocity, sec - ^ 

Subscripts 

c property of the chord 

s property of the shock 

® free-stream condition 

2. INTRODUCTION 

The ability to predict the far-field acoustic signature of a 
high-speed rotor blade has met with limited success in the transonic 
regime. To date, researchers have used the Ffowcs-Williams and Hawkings 
Equation to predict rotor noise in the far field. This method requires 
detailed data for velocity and pressure in the region surrounding the 
blade tip. This information was obtained from finite-difference poten- 
tial solutions of the flow field around a blade rotating at transonic 
tip speeds. This effort produced results which correlated well compared 
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with experimental rotor data taken in high-speed hover. These results 
were valid up to some limiting tip Mach number beyond which the results 
do not correlate well (refs. 1-3). 

When viewed from a blade-fixed coordinate system (fig. 1), the 
velocity of the free stream increases linearly with radius. At some 
radius (defining the "sonic circle"), the free stream becomes sonic. 

The region outboard of the sonic circle is supersonic relative to the 
blade, while the region inboard is subsonic. At high tip Mach numbers, 
a region of supersonic flow, terminated by a shock, will form adjacent 
to and extend away from the blade surface m three dimensions. As the 
tip Mach number increases, the region of supersonic flow will grow 
larger until it finally connects with the outer supersonic region. When 
this occurs, a hyperbolic region will extend from the blade surface to 
the far field. Because of the highly radiative nature of this region, 
the shock will also escape to the far field, dramatically increasing the 
acoustic signature of the rotor. This phenomenon is known as "delocali- 
zation" (ref. 4), and it occurs at subsonic tip Mach numbers. The Mach 
number at which delocalization occurs appears to be the limiting Mach 
number at which theoretical prediction of the far-field acoustic signa- 
ture no longer correlates well with experimental data. It is believed 
that the onset of the characteristic, impulsive "popping" noise of a 
rotor in high-speed flight is caused by this delocalization. 

Delocalization has been investigated using a more basic 2-D model 
(ref. 5). A disturbance in the form of a circular-arc bump spanning the 
surface of a rotating cylinder is used as a computational model. The 
2-D rotational-disturbance model maintains the same mechanism for propa- 
gation of the shock to the far field as that of the 3-D case. The 
relationship of the model to the 3-D rotor case can be seen in figure 2. 
The 2-D computational model is less complex to program and requires much 
less computer time to run. 

As previously stated, velocity and pressure data in the region 
surrounding the blade must be obtained to determine the far-field signa- 
ture. This means that the accuracy of the far-field solution depends 
directly on the accuracy of near-field data. The highly nonlinear 
nature of the flow field at the tip region of the blade requires the use 
of finite-difference methods to obtain the required data. This flow 
field contains weak shocks and the accurate representation of these 
shocks presents a significant problem. Finite-difference formulations 
tend to smear the shock over several grid points. This effect is more 
pronounced with increasing distance from the disturbance. If the shock 
is oblique to the grid, it can become smeared over many points. These 
problems make obtaining an accurate representation of the strongest 
source of far-field noise extremely difficult, and hence the prospects 
of accurate noise prediction are dismal. 
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3. MATHEMATICAL MODEL 


The 2-D computational model is formulated using potential theory. 
This model minimizes computational requirements, while maintaining the 
essential propagation properties. Strictly speaking, the potential 
model does not conserve momentum in transonic problems. However, the 
shocks involved in this problem are typically so weak that the errors so 
caused are negligible — especially in the far field. Furthermore, the 
computational efficiency of a potential model is far greater than that 
of an Euler equation model. The small-disturbance approximation may be 
used to further simplify the model, even though it may be limited at 
some radial distance from the surface because it assumed small shock 
angles. 


The governing equation used as the computational model is 
obtained by casting the potential equation in a reference frame that 
rotates with the circular-arc bump (fig. 3). When the classical small- 
disturbance approximation is invoked, only the lowest-order, nonlinear 
term is retained and the final equation may be written as follows 
(ref. 5) 


where 



(Y + UMjj 





0 


(la) 



co 





The nonlinear sonic circle will be defined at the value of y where the 
coefficient of 4> xx changes sign. 


An interesting feature of this rotational transonic small- 
disturbance equation is the 1/y^ dependence of the nonlinear term. As 
the equation is used at points further away from the surface, this term 
approaches zero and the equation becomes linear. The next higher-order 
term which could be retained results in a modified version, equa- 
tion (1b). This modified equation remains nonlinear in the far field. 
However, tests of this equation show no significant difference from the 
results of equation (la). 



(y 


1)m r 




2M^<J> 

R y xy 


= 0 


(1b) 


4. SHOCK FITTING VS. SHOCK CAPTURING 


The usual method of solving equation (la) is a method known as 
shock "capturing." This method is very desirable because there is no 
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explicit allowance made for the existence of a shock m the flow field. 
Rather, local-differencing schemes are employed to ensure stability in 
both elliptic and hyperbolic regions. During this process, shocks 
appear as high gradients, but not necessarily discontinuous parts of the 
solution. The result is that the shock is smeared over several grid 
points. The smearing occurs because the finite-difference schemes are 
based on Taylor-senes expansions and can be valid only for continuous 
flow regions. Differencing across a shock will result in a locally 
incorrect approximation of the differential equation. For supersonic/ 
subsonic shocks, the smearing will occur over three grid points, while 
for supersonic/supersonic shocks, the smearing may occur over as many as 
ten grid points. The shocks on a typical transonic fixed-wing surface 
are sufficiently strong so that the resolution over three grid points 
exceeds what can be experimentally measured. Therefore, smearing in 
this solution is not objectionable. However, for the problem of high- 
speed rotor noise, the radiative region of interest includes weak 
supersonic/supersonic shocks. This problem is further aggravated by the 
widening of the cylindrical grid as radius increases. The only way to 
improve the resolution of the shock is to add more points to the grid. 
This considerably increases the computation time and the computer 
storage requirements. 

The transonic, small-disturbance equation can be cast in either 
conservative or nonconservative form. Equation (1) is the nonconserva- 
tive form. This form is easier to implement, but mass is not conserved 
across the shock. The reason for this is the overlapping of differenc- 
ing schemes when the shock forms the boundary between a supersonic 
region and a subsonic region. Stability of the numerical method 
requires that a backward- or upwind-differencing scheme be used in a 
hyperbolic flow region, while a central-differencing scheme is used in 
an elliptic region. To conserve mass across the shock, care must be 
taken to write and solve equation (la) in the conservative form as 
follows 
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( 2 ) 


An alternative to the shock capturing method is to explicitly 
impose the shock-jump relation allowed by the transonic small- 
disturbance equation in the solution process. This procedure is known 
as shock "fitting." Such a method, formulated for the potential equa- 
tion by Hafez and Murman (ref. 7) is used m this study. This method 
relies on the a priori knowledge of the presence of a shock in the flow 
field, as well as its approximate location. The implementation of this 
method is outlined in the following steps: 

1. The presence of a shock, as well as its approximate location, 
is determined from a captured solution, either conservative or noncon- 
servative. The location of the shock is determined based on the loca- 
tion of the maximum positive gradient of the pressure coefficient on 
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each grid line of constant y. With the location known, the slope of 
the shock is also determined. From this initial solution, it is also 
possible to determine whether or not the shock has delocalized. 

2. Solve the governing equation using the same capturing method 
up to the shock location. At the first grid point after the shock 
(shock point), solve for the velocity potential explicitly by forcing 
the potential to be continuous across the shock. Also assume a linear 
potential distribution from the shock location to the first computed 
point after the shock point. The value at this grid point will act as 
an internal boundary condition for solution of the equation downstream 
of the shock, which is computed using the capture method. 

3. Once the flow-field solution is obtained, adjust the shock 
location using an unsteady geometric version of the shock-jump relation. 
This will ensure, before the next iteration, that the shock location 
agrees with the latest flow-field solution. Eventually, the shock will 
no longer have to be adjusted and a converged solution will be reached. 

4. Impose continuity of the velocity potential across the shock. 
The movement of the shock during step 3 may result in a downstream grid 
point becoming an upstream grid point or vice versa. If this is the 
case, the value must be corrected by extrapolation to ensure a continu- 
ous potential distribution. 

5. Solve the governing equation by iteration until a steady- 
state shock location is determined. 

Shock fitting conserves mass across the shock because the shock 
is moved so that the flow solution satisfies the shock jump relation. 

As will be seen, this relation is derived from the conservation form of 
the equation. An important rule in implementing the shock-fitting 
method is not to difference across the shock from either side. 

4. 1 Shock-Jump Relation 

The shock-jump relation allowed by the transonic small- 
disturbance equation is the key element in the application of the shock- 
fitting algorithm. Equation (2), which is written in conservative form, 
is of the form 


V - § = 0 

where 

§ = 3 + S 
x y 
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Making use of the shock-tangency condition, and applying the divergence 
theorem, the resultant jump condition is 
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where < > represents the average across the shock. 

When implementing the shock-jump relation, it is necessary to 
give special treatment to certain points in the vicinity of the shock 
location (ref. 6). Moving in the streamwise direction, the last point 
prior to the shock location could look like that shown in figure 4. As 
can be seen, the computational molecule for this point contains the 
point A, which lies downstream of the shock. So as not to violate the 
"rule of forbidden signals," it is necessary to use a value for the 
potential at point A which is extrapolated from the upstream side of 
the shock. Then, the point ahead of the shock may be computed in the 
usual upwind manner. If the shock is supersonic/supersonic, then the 
point after the shock point, B, must be differenced downstream as shown 
in figure 5. Once again, care must be taken not to difference across 
the shock. Since the shock location is known, and 4> must be continu- 
ous across the shock, the value of <j> at the shock is extrapolated from 
upstream data. Just downstream of the shock, <t> x has already been 
computed and serves as a Neuman boundary condition. No special treat- 
ment need be paid to this point if the shock is supersonic/subsonic. 

4.2 Shock Movement 


The requirement to update the shock location requires a form 
of the jump relation in which the shock motion is a function of the 
upstream and downstream flow conditions (ref. 7). Such a relation can 
be obtained from the unsteady small-disturbance equation. 
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The divergence theorem is applied and using the relation 
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the time dependency of 4) can be replaced by a spacial differencing, 
where [[ ]] represents the difference across the shock. This time 
factoring out <j> ultimately results in the following geometric 
relation 


2M 


2 ax 
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This is of the form 


( 7 ) 
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+ B 



C 


(8) 


It can be solved using a variety of methods, but m this instance, an 
implicit trapezoidal scheme is used to avoid any restriction on the size 
of the time step. 




C At B At n+l"l /. B At \ 
A + A Ay x k-1J ‘ V + A Ay/ 


(9) 


Even so, there is an upper limit set not by the stability of the equa- 
tion, but by the amount of shock movement allowed by each iteration. A 
large movement of the shock causes instability m the next iteration of 
the solution. Since the time step is not real, At can be set to any 
value desired by the user. Here, At = 0.5 is used. Once the new shock 
location is found, it is necessary to use a data smoothing technique, 
such as a least-squares fitting, to ensure a continuous shock slope. A 
key feature of equation (8) is that the shock moves based on the right- 
hand side of the equation, which is dependent on the solution to the 
flow field upstream and downstream of the shock. From equation (4), the 
right-hand side of equation (8) defines the inverse of the shock slope 
squared. As the iterations are continued, the slope defined by the 
right-hand side of equation (8) should become equal to the slope on the 
left-hand side. When this has occurred, the solution is said to be 
converged. 

The movement of the shock is the most difficult portion of the 
shock-fitting method to implement. The slope of the shock must be kept 
continuous in some manner. The oblique portion of the shock beyond the 
sonic circle exhibits extreme sensitivity to shock angle. In one case, 
the difference in shock angle of 0.4 of a degree was enough to totally 
destroy the solution. This was a result of the shock- jump relation 
yielding a shock slope approaching the slope of the linear characteris- 
tic as y increases. When this happens, the coefficient of <4> x > in 
equation (4) becomes very small and any small deviation in the shock 
slope from the characteristic value results in large values for <4> x >. 
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5. RESULTS 


To ensure that the shock-fitting method was implemented cor- 
rectly, a step-by-step progression of test cases was accomplished. The 
first case consisted of a wedge in supersonic rectilinear flow with a 
free-stream Mach number of 1.15. This case confirmed the logic neces- 
sary to fit a supersonic/supersonic shock. The result, shown in fig- 
ure 6, was an attached shock at the leading edge. A glance at the 
pressure distribution shows a remarkable clarity of the shock as com- 
pared to the captured case. In the second case, the result is a 
detached bow shock for a circular-arc bump translating at a Mach number 
of 1.15. Figure 7 shows the shock location displacement from the cap- 
tured solution, while figure 8, again, demonstrates the resolution of 
the shock. Note that the solution also shows a shock at the trailing 
edge. The fitting method was not used on this shock. This case is 
important because it demonstrated the shock structure present for the 
rotating, delocalized case which has a supersonic/subsonic portion and a 
supersonic/supersonic portion. 

Once confidence m the implementation of the shock fitting method 
was gained from the preceding test cases, the method was applied to the 
model rotational problem. The rotating circular-arc bump is depicted by 
the following conditions: 1) thickness is 6$, 2) radius/chord ratio is 

10, and 3) rotational Mach number is 0.9. Both Alternating Direction 
Implicit (ADI) and Successive Line Over Relaxation (SLOR) schemes were 
tried and worked successfully, with ADI having the expected convergence 
rate advantage. Figure 9 shows the comparison of the captured shock 
location with that of the fitted solution. Since the captured solution 
is formulated using the nonconservative form, the shock is slightly 
forward of the fitted solution which conserves mass across the shock. 
This is to be expected within the sonic circle. Beyond the sonic 
circle, there is no conservation question since there is no need to use 
a switched-differencing scheme. The shock moves forward m the course 
of correcting the extreme smearing which is exhibited by the captured 
solution. Figures 10-12 show the comparison of captured and fitted 
pressure data at different radii from the surface of the cylinder. In 
figure 10, the results show very little difference in amplitude or in 
the discontinuous nature of the shock. As the distance from the distur- 
bance increases (figs. 11 and 12), the solutions become much different. 
The captured solution shows marked smearing, while the fitted solution 
not only shows a sharp discontinuity, but also marked increase in the 
amplitude of the pressure. 

Figure 13 gives an overview of the flow field and compares the 
captured and fitted solutions by using isomach lines. Here, it is clear 
that the shock maintains a sharp discontinuity throughout the flow 
field. 
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6. CONCLUSIONS 


The 2-D computational model using the TSD equation models the 
delocalization phenomenon well. It is simple in all computational 
aspects, yet still provides the mechanism which allows shocks to propa- 
gate to the far field from a subsonically rotating disturbance. Apply- 
ing the shock-fitting algorithm to the solution process changes the 
nature of the resulting flow field significantly. The discontinuous 
nature of the shock is maintained throughout the field and well beyond 
the sonic circle. In the process of sharpening this discontinuity, the 
amplitude of the shock increases. Both of these factors are important 
in predicting the acoustic signature in the far field, since the shock 
is the primary near-field noise source. 

The implementation of the shock-fitting scheme for the 3-D rotor 
case is being pursued. Recently measured, experimental rotor data will 
be used to verify the computed near-field acoustic signature. This will 
validate the accuracy of the computed, near-field data necessary to 
predict the far-field signature. 

The shock-fitting method has some minor drawbacks. The complex- 
ity of programming is increased because of the logic necessary to accom- 
modate the shock points. A stretched grid may still limit the resolu- 
tion of a shock far away from the disturbance; however, the solution 
will always show a considerably sharper shock than the captured method. 
The TSD equation may also be of limited use in the far field, but it may 
be adequate for determining the desired flow properties in the near 
field. 


When applied to the 3-D case, the solution should provide insight 
into the theoretical prediction of the far-field acoustic signature for 
a rotor operating at tip Mach numbers beyond the delocalization Mach 
number . 
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Fig. 1. Blade-fixed coordinate system from reference 4. 
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Comparison of computational model to actual rotor. 
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Coordinate transformation of computational model. 



Fig. 4. Differencing upstream of a shock. 
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Fig. 5. Differencing downstream of a supersonic/supersonic shock. 
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Fig. 6. Comparison of fitted and captured solutions for a wedge at Mach 
number = 1.15. 



Fig. 7. Comparison of shock locations for a 6$ circular arc translating 
at Mach number = 1.15. 
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TRANSLATING CIRCULAR ARC 
VERTICAL LOC =3 03 CHORDS 
MACH NO = 1 15 
THICKNESS = 6% 



Fig. 8. Comparison of fitted and captured solutions for a 6 % circular 
arc translating at Mach number = 1.15. Solution is 3.03 chords above 
the surface. (Rear shock is not fitted.) 



Fig. 9. Comparison of shock locations for a 6 % circular arc rotating at 
Mach number = 0.9. 
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CHORDS FROM CENTER 


Fig. 10. Comparison of fitted and captured solutions for a 6% circular 
arc rotating at Mach number = 0.9; y = 0. 



Fig. 11. Comparison of fitted and captured solutions for a 6$ circular 
arc rotating at Mach number = 0.9; y = 3.24 chords. 
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VERTICAL LOC - 8 85 CHORDS 
ROTATIONAL MACH NO = 0 9 



Fig. 12. Comparison of fitted and captured solutions for a 6% circular 
arc rotating at Mach number = 0.9; y = 8.85 chords. 
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MACH = 1 00 


Fig. 13. Isomach contours comparing a captured solution and a fitted 
solution for a 6 % circular arc rotating at Mach number = 0.9. 
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